function dF = ffunODE(V,F, r, gamma ,Lambda, phi,q,aBar,delta)
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Solve HJB Equation %%%%%%%%%%%%
%Calculate second derivative F''(V), given F(V), F'(V), and V
% Applies for general case with maturity, i.e., delta>0
% Baseline is obtained for delta=0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        dF=zeros(2,1);
        dF(1)=F(2);
        
        F1=F(1);
        
        %%%%%% Calculate Monitoring Effort according to FOC
        a=(F1-dF(1)*(V+phi)-phi*(gamma-r))/phi;

        %%%%%%%%% Impose lower and upper boundary
        if (a<0)
            
            a=0;
          
            
        end
        
        if (a>aBar)
            
            a=aBar;
           
            
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%% Impose that W=a\phi<F(V) %%%%%%%%%%%%%
        %%%%% If effort from FOC exceed, F(V)/\phi, set a=F(V)/phi
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        if (phi*a>F1)
            
            %Case F(V)=W(V)

            a=F1/phi;
            lambda=Lambda-q-a;
            W=phi*a;
            dV=(gamma+lambda+delta)*V-W;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%% Set Second Derivative %%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            dF(2)=(dF(1)^2-dF(1)*F1/phi+dF(1)^2*V/phi-delta*dF(1))/(min(dV,-0.0001));

            %Impose lower boundary to avoid degenerate solution; this
            %constraint does not bind in optimum and is for numerical
            %purposes only
            dF(2)=min(dF(2),-0.0001);

        else

            %Case F(V)>W(V)

            lambda=Lambda-q-a;
            W=phi*a;
            dV=(gamma+lambda+delta)*V-W;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%% Set Second Derivative %%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            dF(2)=-(gamma+delta-r)*dF(1)/(min(dV,-0.0001));

            %Impose lower boundary to avoid degenerate solution; this
            %constraint does not bind in optimum and is for numerical
            %purposes only

            dF(2)=min(dF(2),-0.0001);
            
        end

end